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ALGORITHMS  FOR  INTERPOLATING  AND  APPROXIMATING 
RADIANCE  PROFILES 

Sieven  J.  Leon 


1  Introduction 

In  recent  years  there  has  been  a  considerable  amount  of  work  by  J.  I.  F.  King  and  others 
developing  methods  for  remote  temperature  sensing  based  on  applying  transform  theory 
to  the  radiative  tranfer  equation.  King,  Hohlfeld,  and  Kilian,  [9]  have  shown  that  using 
differential  inversion  it  is  possible  to  successfully  determine  temperature  profiles  based  on 
measurements  of  the  upwelling  intensities  in  the  atmosphere.  A  second  transform  method 
proposed  by  Jean  I.  F.  King  is  based  on  optical  measure  theory.  This  technique  has  been 
studied  by  King  and  Leon,  [11,12].  The  optical  measure  theory  method  has  the  advantage 
that  i'  djes  not  require  any  computations  of  numerical  derivatives.  Instead  the  radiance 
profil  '.s  approximated  by  a  rational  function.  The  Planck  intensity  is  then  determined 
as  t.ij  inverse  transform  of  the  rational  function.  In  this  paper  we  present  techniques 
for  aoproximating  the  radiance  profile  based  on  rational  interpolation  and  nonlinear  least 
scuaics  data  fitting. 

In  the  final  section  of  the  paper  we  discuss  general  matrix  methods  for  solving  integral 
equations  of  the  first  kind.  In  particular,  using  regularization  techniques,  the  inversion 
problem  can  be  solved  as  a  linear  system  with  quadratic  constraints.  Numerical  algorithms 
have  been  developed  for  solving  these  constrained  problems.  However,  these  algorithms 
proved  unsuccessful  when  applied  to  data  sets  obtained  from  the  NOAA  TIROS  Operational 
Vertical  Sounder  (TOYS).  A  singular  value  analysis  is  given  which  shows  why  the  matrix 
techniques  fail. 


2  Rational  Approximation  of  Radiance  Profiles 


Using  radiative  transfer  theory  one  can  relate  the  Planck  intensity  to  the  upwelling  intensity 
of  the  atmosphere.  This  relationship  can  be  expressed  in  terms  of  an  integral  equation  of 
the  first  kind.  Specifically  the  relationship  is  given  by 


/?(p)  =  /  D(p)W(p/p)dp/p  (1) 

Jo 

where  W{p/p)  is  a  kernel  weight  function  that  peaks  at  p  =  p  and  R(p)  denotes  the  radiance 
of  the  wavelength  channel  whose  weight  function  peaks  at  p  =  p. 

If  we  set 

W(r)  =  ze"*  and  s  =  - 

P 


□o 


then  the  integral  equation  (2)  can  be  expressed  as  a  Laplace  transform  and  consequently  we 
can  solve  for  the  Planck  function 

We  can  think  of  B{p)  and  R{p)  as  transform  pairs.  Once  a  representation  for  R  has  been 
decided  upon,  then  B  can  be  determined  analytically  as  an  inverse  Laplace  transform.  In 
Section  4  we  will  consider  other  choices  for  the  weight  function  VV{c).  Ontical  measure 
theory  permits  the  inversion  of  non-Laplacian  exponential  kernels. 

The  ^-function  inversion  theory  of  Chandrasekhar  [2]  suggests  that  R  should  be  repre¬ 
sented  as  a  rational  function.  If  we  set 


then 


and 


R{p)  —  dj  -b  dop-i - 1-  -b 

i=\  ^ 

s  s  ^  ®  +  <^> 


(2) 


B{p)  =  di  +  d'ip  -b  •  ■  •  + 

i  =  l 


-C.p 


(3) 


The  coefficients  in  equation  (2)  can  be  determined  from  the  data  by  rational  interpolation 
or  by  nonlinear  least  squares.  If  one  sets 


R,  =  R{pi)  »=I . m 


and  expresses  fZ  as  a  quotient 


F  -  an+2P'  ‘  -  •  •  •  -  Om 
then  R  will  interpolate  (p,-,  Ri)  if  and  only  if 

“iPi”  +  - b  Oi»+i  +  ar>+2p/”'  -b - b  =  Pi'Ri 


(1) 


i= 

The  coefficients  a*  are  determined  by  solving  the  linear  system  (4).  The  residue  form  of  R 
can  be  computed  as  an  inverse  convolution  and  from  this  one  can  determine  the  coefficients 
for  equations  (2)  and  (3).  In  practice  these  equations  will  only  have  a  small  number  of 
components.  Generally,  it  is  possible  to  get  a  good  fit  with  j  =  1  or  2  and  /  =  3  —  y.  If  one 
takes  i  =  1  and  /  =  2,  it  is  possible  that  for  some  data  sets  the  interpolating  function  will 
have  a  positive  pole  (cj  <  0),  even  though  this  is  impossible  for  an  actual  radiance  profile. 
When  this  happens  the  weight  Wj  corresponding  to  the  pole  will  be  much  smaller  than  the 
weight  corresponding  to  the  other  rational  component.  Thus  if  one  sets  Wi  =  0,  then  the 
resulting  function  will  give  a  good  approximation  to  the  data  points. 
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The  rational  interpolating  function  can  be  taken  as  starting  approximation  for  an  itera¬ 
tively  computed  nonlinear  least  squares  fit  to  the  data  of  the  form 


R(p)  =  xip-f  xj-f 


X3 

1  +  Xi^p  1  +  xs^p 


(5) 


Initially  the  Xi  coefficients  are  determined  by  the  coefficients  of  the  rational  interpolating 
function.  If  the  interpolating  function  has  no  positive  poles,  then  set 


xi=ai,  xi  =  ai,  X3  =  u'i,  X4  =  w^,  xs  =  y/cl,  x^  =  y/c^ 


If  the  first  pole  of  the  interpolating  function  is  positive,  set  X3  =  0  and  X5  =  1.  Similarly  if 
the  second  pole  is  positive  set  x^  =  0  and  xg  =  1. 

The  choice  of  a  rational  function  to  represent  the  radiance  data  is  motivated  by  the 
physics  of  the  atmosphere.  In  practice  the  nonlinear  least  squares  approximation  gives  a 
very  good  fit  to  the  radiance  data.  The  coefficients  of  the  rational  function  can  be  used  to 
determine  a  Planck  function  of  the  form  (3).  ! 

■  j 

3  Generalized  Exponential  Inversion 

In  the  generalized  exponential  inversion  proposed  by  Jean  I.  F.  King  it  is  a.ssumed  that 

H'(..)=VFi(.-)  =  7i.'e.xp*(c)  ! 

I 

where  i 

= _ ! _  i 

and  I 

e.xpt(i)  =  e\p[-:'‘/k) 

King  has  shown  that  if 

H(p)=  r 

Jo  P 

and  R{p)  is  of  the  form  (2).  then  the  Planck  intensity  D(p)  =  Bk{p)  will  be  of  the  form 


d 

Bkip)  =  di  -t-  dnp  +  •  •  •  +  -^p’  +  ^  expt(rip) 

i=l 

Note  that  equations  (3)  and  (6)  are  the  same  in  the  case  fc  =  1. 


{'5) 


4  Test  Results 

The  interpolation  and  nonlinear  least  squares  algorithms  have  been  tested  extensively  on 
simulated  data  and  on  the  TOYS  data.  The  algorithms  are  numerically  stable.  Figures 
1  through  5  were  generated  from  a  typical  data  set  of  the  TOYS  data.  In  Figure  1  the 
radiance  has  been  derived  by  a  nonlinear  least  squares  fit  of  the  form  given  in  equation  (-5). 
The  radiance  values  are  plotted  as  a  function  of  the  pressure  in  millibars.  Figure  2  shows 
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a  semilog  plot  of  the  temperature  profile  for  the  same  data  set.  The  pressures  are  given 
in  millibars  on  the  vertical  axis.  Figure  3  shows  the  temperature  profile  for  that  data  set 
derived  by  generalized  exponential  inversion  with  k  =  0.8.  Figure  4  shows  the  temperature 
profile  derived  by  generalized  exponential  inversion  with  k  =  1.2.  Figure  5  shows  all  three 
Planck  plotted  as  functions  of  p  on  the  same  axis  system. 


5  Matrix  Methods 

Equation  (1)  can  be  regarded  as  a  Fredholm  integral  equation  of  the  first  kind.  By  this  we 
mean  an  equation  of  the  form 

J  f{y)K(x,y)dy  =  g{x)  (7) 

where  the  functions  g{x)  and  K{x,y)  are  known  and  the  solution  f(y)  must  be  determined. 
The  difficulty  v/ith  this  type  of  equation  is  that  the  solution  does  not  depend  continuously 
on  the  data.  U  is  possible  to  have  function  H{y)  that  oscillates  wildly  on  [a,  6]  such  that 

H(y)I<{x,  y)dy  =  f(i)  «  0 


Thus  if  F  is  a  solution  to  (7),  then  Fi  =  F  +  H  is  a  solution  to 

J  /(y)A'{Al/)dy  =  <?i(ir) 

where  <7i(a:)  =  g[x)  +  e{x).  The  functions  F  and  Fj  will  differ  greatly  even  though  g 
and  Qi  are  close.  In  practice  the  known  function  g{x)  is  represented  by  a  discrete  data 
set  which  involves  some  experimental  error.  In  general  because  the  problem  is  ill-posed, 
the  numerically  computed  solution  is  not  likely  to  be  close  to  the  true  solution.  This  is 
particularly  true  when  the  function  g{x)  is  sampled  at  only  a  small  number  of  data  points. 

In  this  section  we  will  discuss  general  matrix  methods  involving  regularization  conditions 
which  can  be  used  for  solving  integral  equations  of  the  first  kind.  We  will  examine  how  well 
these  methods  can  be  applied  to  the  transfer  equation. 

Equation  (7)  can  be  discretized  using  an  appropriate  quadrature  formula  to  give  an 
equation  of  the  form 

n 

;=i 

Setting 


Oij  =  WjK(xi,yj) 


-j=/(yj)  l>i  =  g(xi)  f=l,...,m  j=l,...,n 


the  discretized  equation  can  then  be  represented  as  a  linear  system 


Az  =  h 
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If  A  has  singular  value  decomposition  then  the  solution  is  given  by 

2  =  vl+b  = 

Because  the  problem  is  ill-posed,  the  computed  solution  z  may  vary  greatly  from  the  true 
solution.  A  standard  technique  is  to  require  that  the  solution  /(y)  satisfy  a  smoothness 
constraint  in  order  to  eliminate  wildly  oscillating  components.  The  regularization  condition 
one  usucdly  imposes  is  that  the  norm  of  one  or  more  of  the  derivatives  of  /  be  bounded. 
When  this  constraint  is  discretized  one  ends  with  a  constrained  least  squares  problem  of  the 
form 

i.Iinimize 

||Az  -  bli2 

subject  to 

lICzi|j<<i  (CeR'”'") 

This  constrained  problem  can  be  solved  using  the  method  of  W.  Gander  [5].  First  we 
note  that  if  .d+b  is  not  feasible,  then  the  solution  z  must  satisfy  the  constraint  l|Cz||2  =  d. 
We  are  then  left  with  the  equality  constrained  problem  which  ca.i  be  solved  using  Lagrange 
multipliers.  This  leads  to  the  generalized  normal  equations 

(A^A  +  AC^C)z  =  A^b 

These  equations  can  be  solved  using  the  generalized  singular  value  decomposition, 

f/^AA' =  diag(oi,.,.,an)  =  f?! 

=  diag(/?i . /J,)  =  D2 

where 

q  =  min(p,  n)  U'^U  =  V^V  =  Ip 

and  X  is  a  nonsingular  n  x  n  matrix,  if  we  set  w  =  A’“'z  and  c  =  U^h,  then 

||Az-b|l2  =  llD,w-c||2 

and  the  generalized  normal  equations  reduce  to 

{DjDi  +\DXD2)yf  =  DJc 

We  can  solve  for  w  as  a  function  of  A. 

«*^AA)  =  -j  ^  j=I,...,n 

-I-  A/i/ 

The  secular  equation 

v'(A)  =  ||D2w(A)||2  =  d^- 

will  have  only  one  positive  re,al  solution  A*  and  this  is  the  solution  which  minimizes 
llDiw(A)  —  c||2.  (See  Gander  f.5]).  The  solution  to  the  constrained  least  squares  problem  is 
then  given  by 


z  =  .Yw(A*) 


In  the  case  that  the  system  i4z  =  b  is  underdetermined  (m  <  n)  one  may  consider 
reversing  the  conditions.  This  leads  to  the  second  regularization  problem: 

Minimize 


subject  to  the  condition 


||Xz  -  b||j  =  min 


This  second  problem  is  a  much  easier  problem  to  solve.  If  A  has  singular  value  decomposition 
f/EV’’  and  .  . 


[o  0) 


where  Ei  is  nonsingular,  then  set 


We  can  minimize  \\Az  —  b((2  by  setting 

yi  =  Ef'ci 

The  vector  yj  must  then  be  chosen  so  as  to  minimize  HCzllj.  To  do  this  set 


It  follows  that 


E  =  CV  =  (£]  Ej)  and  e  =  -£iyj 


r»ll2  =  \\Ey\\2 

=  ll^’zyz  “  ®ll2 


Thus  l|Cz|l2  can  be  minimized  by  setting  y2  =  E^e.  The  solution  to  the  second  regular¬ 
ization  problem  is  then  given  by 

z  =  V'y 

These  regularization  iiKthods  can  be  applied  to  the  transfer  equation.  Generally,  because 
the  only  useful  measurements  are  those  taken  at  wavelengths  near  15  /rm  ,  one  is  usually 
limited  to  approximately  6  data  points.  Thus  if  equation  (1)  is  discretized  using  16  point 
Gauss-Laguerre  quadrature,  then  the  resulting  coefficient  matrix  A  will  have  dimensions 
6  X  16.  In  practice,  most  of  the  singular  values  of  /I  will  be  quite  small.  For  a  typical  matrix 
determined  by  a  sample  of  data  from  the  NOAA  TIROS  Operational  Vertical  Sounder 
the  computed  singular  values  (rounded  to  2  digits)  were:  0.28  x  10°,  0.45  x  10“*,  0.60  x 
10“*,  0.22  X  10“®,  0.30  X  10“*,  0.78  x  10“**.  Because  the  numerical  rank  of  A  is  so  low 
relative  to  the  precision  of  the  data,  the  matrix  methods  discussed  earlier  do  not  work.  If 
the  first  matrix  method  is  used  then,  for  any  reasonable  constraints,  the  computed  residuals 
turn  out  to  be  unacceptably  large.  On  the  other  hand,  if  the  second  method  is  used,  it 
is  not  possible  to  obtain  a  reasoi:abie  constraint  at  all.  There  just  simply  is  not  enough 
information  for  the  regularization  methods  to  work. 
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6  Conclusions 


Although  transfer  theory  can  be  used  to  relate  the  Planck  Intensity  to  the  upwelling  intensity 
in  the  atmosphere,  the  relation  is  expressed  in  the  foim  of  an  integral  equation  of  the 
first  kind.  Such  equations  are  ill-posed  and  consequently  do  not  have  unique  numerical 
solutions.  Adding  regularization  conditions  does  not  solve  the  problem.  For  data  sets 
like  those  obtained  from  the  TOYS  data,  there  is  not  enough  information  to  determine 
the  Planck  function  using  matrix  methods  with  appropriate  smoothness  constraints.  More 
assumptions  relative  to  tlie  physics  of  the  atmosphere  must  be  added.  In  this  regard  the 
transform  methods  based  on  optical  measure  theory  seem  to  work  well.  The  key  step  is  the 
representation  of  the  radiance  profile  by  a  rational  function.  This  can  be  accomplished  in  a 
numerically  stable  manner  using  nonlinear  least  squares  fitting  algorithmns. 
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